########################################################################
#Replication data for
  #A replication of "Exploring and explaining contracting out: Patterns among the American States"
#By: Jing Qian, Jiahuan Lu, Jianzhi Zhao
#Date: 09/06/2022
########################################################################


########################################################
#TABLE 2
########################################################
#---------------
#(1) Setup
#---------------
rm(list = ls())

library(tidyverse)

#---------------
#(2) Load data
#---------------
#Raw ASAP survey data
load("replication_data_table2.RData")

#Only use 1998, 2004, 2008
asap.use = asap.raw %>%
  subset(year %in% c(1998, 2004, 2008))

#Only include those use contracting out
asap.contract = asap.use %>%
  subset(l_1a == 1)

#---------------
#(3) Calculate ratios
#---------------
#Engagement in contracting
t2.1 = asap.use$l_1a %>%
  tapply(., asap.use$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,2]
  
#Budget: 5% or less
t2.2 = asap.contract$l_1b %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,1]

#Budget: 40% or more
t2.3 = asap.contract$l_1b %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,4]

#Contract with: other governments
t2.4 = asap.contract$l_2a %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,2]

#Contract with: nonprofit organizations
t2.5 = asap.contract$l_2b %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,2]

#Contract with: for-profit businesses
t2.6 = asap.contract$l_2c %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  .[,2]

#Effect on quality
t2.7 = asap.contract$l_3a %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  lapply(., function(x) x[1:3]) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  t()

#Effect on cost
t2.8 = asap.contract$l_3b %>%
  tapply(., asap.contract$year, table) %>%
  lapply(., prop.table) %>%
  lapply(., function(x) x[1:3]) %>%
  do.call(rbind, .) %>%
  round(., digits = 3) %>%
  (function(x) x * 100) %>%
  t()

#---------------
#(4) Combine
#---------------
#Combine
t2.combined = rbind(t2.1, 
                    t2.2,
                    t2.3,
                    t2.4,
                    t2.5, 
                    t2.6,
                    t2.7,
                    t2.8) %>%
  data.frame()

#Add rownames
t2.combined = cbind(c("Engagement in contracting",
                      "Contracting as % of agency budget: 5% or less",
                      "Contracting as % of agency budget: 40% or more",
                      "Contract with: other governments",
                      "Contract with: nonprofit organizations",
                      "Contract with: for-profit businesses",
                      "Effect on quality: no effect",
                      "Effect on quality: increased",
                      "Effect on quality: decreased",
                      "Effect on cost: no effect",
                      "Effect on cost: increased",
                      "Effect on cost: decreased"), 
                    t2.combined)

colnames(t2.combined) = c("", "1998", "2004", "2008")
rownames(t2.combined) = 1:nrow(t2.combined)

#Display
print(t2.combined)

#---------------
#(5) Save
#---------------
writexl::write_xlsx(t2.combined,
                    path = "result_table2.xlsx")




########################################################
#Table 3
########################################################
#---------------
#(1) Setup
#---------------
rm(list = ls())

library(tidyverse)
library(lmerTest)
library(modelsummary)

#Custom function
make.equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}
#---------------
#(2) Load data
#---------------
load("replication_data_table3.RData")

#---------------
#(3) Column 2: Narrow Replication (1998)
#---------------
#Varlist
dep.var = "contracting"
vars.98.narrow = c("budgexpan",
                   "govcntrl",
                   "agenideo",
                   "fedaction",
                   "budget",
                   "prevcontra_quality_factor",
                   "prevcontra_cost_factor",
                   "experience",
                   "tenure",
                   "education",
                   "sectorpref",
                   "reinvent",
                   "agency",
                   "competition",
                   "costsav",
                   "fiscalrev",
                   "fiscaldemd",
                   "citizpref",
                   "govideo",
                   "pubemplo",
                   "popchang",
                   "popblck",
                   "perfund",
                   "gpphrgrd")


#Make equation
eq.98.narrow = make.equation(dep.var,
                             vars.98.narrow) %>%
  paste0(., " + (1|state)")

#Fit
t3.2 = lmer(formula = eq.98.narrow,
            data = asap.98,
            REML = FALSE)

#---------------
#(4) Column 3: Wide Replication (2004 & 2008)
#---------------
#Varlist
dep.var = "contracting"

vars.0408 = c("budgexpan",
              "govcntrl",
              "agenideo",
              "fedaction",
              "budget",
              "prevcontra_quality_factor",
              "prevcontra_cost_factor",
              "experience",
              "tenure",
              "education",
              "reinvent",
              "agency",
              "competition",
              "costsav",
              "fiscalrev",
              "fiscaldemd",
              "citizpref",
              "govideo",
              "pubemplo",
              "popchang",
              "popblck",
              "gpphrgrd")



#Make equation
eq.0408 = make.equation(dep.var,
                        vars.0408) %>%
  paste0(., " + (1|state)")

#Fit
t3.3 = lmer(formula = eq.0408, 
            data = asap.0408,
            REML = FALSE)
#---------------
#(5) Column 4: Wide Replication with year FE (2004 & 2008)
#---------------
#Make equation
eq.0408.year = make.equation(dep.var,
                             c(vars.0408, "factor(year)")) %>%
  paste0(., " + (1|state)")

#Fit
t3.4 = lmer(formula = eq.0408.year, 
            data = asap.0408,
            REML = FALSE)

#---------------
#(6) Column (5): Wide Replication with two-way FE (2004 & 2008)
#---------------
#Make equation
eq.0408.fe = make.equation(dep.var,
                           c(vars.0408, 
                             "factor(year)",
                             "factor(state)")) 

#Fit
t3.5 = lm(formula = eq.0408.fe, 
          data = asap.0408)


#---------------
#(7) Combine Results
#---------------
t3.combined = modelsummary(models = list(t3.2,
                                         t3.3,
                                         t3.4,
                                         t3.5),
                           coef_omit = "year|state|Intercept",
                           conf_level = 0.90,
                           fmt = 4,
                           stars = c("*" = 0.1,
                                     "**" = 0.05,
                                     "***" = 0.01),
                           statistic = "p.value",
                           output = "data.frame")


#---------------
#(8) Order and Format
#---------------
#Remove extra info
t3.combined = t3.combined %>%
  .[, c(2, 4:7)] %>%
  .[c(1:74, 76),]


#Remove agencies with no sig results
agency.remove = c("agencyStaff-Fiscal",
                  "agencyEducation",
                  "agencyNatural Resources",
                  "agencyCriminal Justice",
                  "agencyRegulatory",
                  "agencyOther")
t3.combined = t3.combined %>%
  subset(! term %in% agency.remove)

rownames(t3.combined) = 1:nrow(t3.combined)

#Add FEs
t3.combined = rbind(t3.combined[1:62,],
                    c("State Random Effects", "Yes", "Yes", "Yes", "No"),
                    c("State Fixed Effects", "No", "No", "No", "Yes"),
                    c("Year Fixed Effects", "No", "No", "Yes", "Yes"),
                    t3.combined[63,])

  
rownames(t3.combined) = 1:nrow(t3.combined)

#Format
t3.combined$term[seq(2, 62, 2)] = ""

colnames(t3.combined) = c("", 
                          "(2) Narrow Replication (1998)",
                          "(3) Wide Replication (2004 & 2008)",
                          "(4) Wide Replication with year FE (2004 & 2008)",
                          "(5) Wide Replication with two-way FE (2004 & 2008)")

#Display
print(t3.combined)

#---------------
#(9) Save
#---------------
writexl::write_xlsx(t3.combined,
                    path = "result_table3.xlsx")














